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We investigate the critical end point of QCD with two flavours of light dynamical quarks at 
finite lattice cutoff a — 1/4T using a Taylor expansion of the baryon number susceptibility. We 
find a strong volume dependence of its radius of convergence. In the large volume limit we obtain 
/is/T ~ 1.1 as the radius of convergence at T/Tc ~ 0.95 where Tc is the cross over temperature at 
zero chemical potential. Since this estimate is a lower bound on the critical end point of QCD, the 
above small value of jj, may place it in the range of observability in energy scans at the RHIC. 
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' The phase diagram of QCD contains two experimentally tunable couplings — the temperature, T, and the baryon 

' chemical potential, /i^. For physical values of the quark masses, it is expected to have a line of first order phase 
transitions in the T-fiB plane. This line rises from the axis to terminate at the critical end point, which is specified 

I by its coordinates, (T^, Hg). There is no critical point on the T axis, but there is a rapid cross over in quantities such 

^ . as the Wilson line or the quark condensate at a temperature Tc. This temperature is therefore identifiable through 

jy-^ ' peaks in the susceptibilities, i.e., the temperature derivatives, of these quantities. In the limit of vanishing quark mass, 

, Tc becomes a critical point and the critical end point becomes a tricritical point S 13 ■ Model estimates of the 
location of the end point vary widely, thus calling for a first principles determination through a lattice computation. 

04 , If /i^ is small enough that it can be reached in present day experiments, then it may be detectable in a variety of 

• ways 0. 

I Direct simulations of QCD at finite chemical potential are not possible due to the fermion-sign problem. However, 
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I. INTRODUCTION 




FIG. 1: The method followed in this work seeks to Taylor expand |3 results obtained at vanishing fi along lines of constant T; 
thereby bracketing the critical region, shown by circles. This region can vary in location (as shown by the dashed line) and size 
(indicated by the circles) as the volume changes, growing smaller with increasing volume. Measurements of these changes can, 
in principle, be used to identify the critical exponents and thereby pin down the universality class at the critical end point. 
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various methods have been developed to continue lattice simulations at computable points in the phase diagram to 
the more interesting and uncomputable regimes. These methods include reweighting and its variants ||6i,i7|], analytic 
continuation of computations at imaginary chemical potential |^ 0| and Taylor series expansions of the free energy 

In this paper we make the first attempt to locate the critical end point on large volumes and small pion masses 
through a Taylor expansion of the quark number susceptibility in the variable /is at several different T. This Taylor 
expansion allows us to extrapolate out in fiB along lines of constant T (see Figure^. Computing the extrapolation 
at several T allows us to bracket the critical region and track its change as the lattice volume is changed. Systematic 
uncertainties in setting the scale of T (see Section Hvj) determine how closely these extrapolation lines can be spaced. 
Typically, these uncertainties decrease with lattice spacing. In the future, as one approaches the continuum limit with 
increased computing power, it should become possible to obtain the critical exponents using this method. 

At this time however, computational resources are not sufRcient to complete this program. In this work we concen- 
trate on bracketing the critical region and investigating gross changes with volume. All finite volume estimates should 
technically be called pseudo-critical values, which tend to the critical values on infinite lattices. On large enough 
lattices the distinction is often immaterial. Our larger lattices are the largest ever used in investigating this problem, 
and the estimates of the radius of convergence that we obtain here are compatible with the infinite volume results, as 
we discuss later. 

Operating within the context of the present understanding of the phase diagram, we scan downwards in temperature 
from T > Tc- The first sign of a phase transition that one meets is then an estimate of the critical end point. 

Our results could be compared with earlier work on the estimation of the end point with either two flavours (up and 
down) of light staggered quarks {Nj = 2) or two light (up and down) and a third heavier (strange) staggered quark 
{Nf = 2-1-1). All such computations have been performed on lattices with temporal extent Nt = 4, i.e., with lattice 
spacing a = 1/4T. All have comparable lattice spacings in physical units, since nip/Tc = 5-5.5. Nevertheless, the 
quark masses and lattice volumes {N^., where is the spatial extent) differ widely (a detailed comparison is given 
later in Table ITTjl . The Nf = 2 computation of 0| used a very high up/down quark mass although the lattice 
size was large in units of the pion's Compton wavelength. Significantly lighter up/down quark masses were used on 
much smaller lattices in the Nf = 2-1-1 study of 0| and the Nf — 2 work of . The only errors reported for 
these estimates are statistical (see Table |n] later). 

In fact, the largest errors are likely to be systematic. The two computations of ^| indicate that quark mass 
effects are large. Finite lattice spacing effects have also been estimated to be as large as the statistical errors Q. 
Moreover, one should expect two kinds of finite volume effects. A strong "small" volume effect must arise from 
distortions of the spectrum of the Dirac operator when the spatial size is as small as 2-3 pion Compton wavelengths. 
These must be removed before any physics is extracted. The more benign "large" volume effects are welcome, since 
it is the analysis of such effects which would eventually verify whether one is indeed studying a critical point. In 
view of this we decided to work with Nf — 2 QCD with intermediate quark masses giving rriT^/mp = 0.31 ± 0.01 at a 
lattice spacing of 1/4T where nip/Tc = 5.4 ± 0.2. We expect our results to be comparable with those of 0] because 
the scales determined at = are comparable. We differ from earlier studies with similiar quark masses in that our 
lattices are significantly larger. From our estimate of the radius of convergence on lattices with spatial sizes in the 
range NsrriT^ « 3-10, we found strong "small" volume effects at the lower end of Ng, and a more stable physical value 
for larger Ng. 

The plan of our paper is as follows — the Taylor expansion is introduced and the quantities to be evaluated on 
the lattice are set out in detail in Section ^] Efficient numerical techniques for performing the trace computations 
are given in Section IIIII along with details of the performance of the algorithms. The determination of simulation 
parameters, the details of the simulations and the numerical results are given in Section IIVI The main results are 
collected and discussed in Section Detailed formulae for the non- linear susceptibilities are given in the apppendix. 
The paper is organized such that a perusal of Section alone would satisfy a reader who is familiar with the context 
and needs to extract our main results. If such a reader finds the language to be unfamiliar, then referring to Section 
m would suffice. 

II. THE TAYLOR EXPANSION 

The partition function for QCD at temperature T and chemical potentials /i/ for each of Nf flavours, can be written 
in the form 

Z{TAh}) = / 2?f/e-^«(^) X{TietMf{mf,T,^if), (1) 
■' f 
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where Sq is the gluon part of the action and M denotes the Dirac operator. We employ the standard Wilson action for 
Sg and staggered fermions to define M. The fermion-sign problem, which prevents direct Monte Carlo simulations 
of QCD at finite chemical potential, is due to the fact that the determinants have arbitrary complex phases for 
non- vanishing /iy. One possible solution to is to recognize that the pressure, 

P{T,M)^-^= (^^yogZiT,M), (2) 

can be expanded in a Taylor series about the point where all the A*/ = 1^. The leading term, independent of /x/, 
has been obtained in previous computations |l3l |. The pressure is a convex function of the intensive thermodynamic 
variables, T and {a*/}- The first derivative of the pressure with respect to fij is the quark number density for flavour 
/. The second derivative is the corresponding quark number susceptibility (QNS) As one approaches the critical 
end point, the diagonal QNS [l5l | diverges as a power-law determined by the critical exponents. The location of the 
divergence gives (T^,^^), and a determination of the critical exponent could verify the arguments that the QCD 
critical end point is in the same universality class as the 3-d Ising model 4] . 

In this paper we deal with Nf — 2 and mainly with the baryon chemical potential fJ-B — 3/x„ — Sfid- The quarks 
have a small but non-vanishing mass, m„ — nid — m. The A^-th order derivatives in the Taylor expansion then 
can be taken n„ times with respect to fiu and nd = N — times with respect to fid- We denote this non-linear 
quark number susceptibility (NLS) as Xn„,na- Explicit operator expressions for the NLS are given in the appendix, 
along with an exposition of the methods and operators, 0„, used to compute them. Flavour symmetry implies that 
Xn^,na — Xna,nu- The Taylor expansion of the chemical potential dependent part of the pressure, 

AP{T, fid) = P{T, fiu, Md) - P{T, 0,0)= V xn„.«, ^ ^ (3) 

can be translated into a joint Taylor expansion in the baryon chemical potential fiB = 3(/iu + Mrf) /2 and the iso- vector 
chemical potential /i/ = ~ fid)/^- This double expansion can be specialized to one in ^bI"^ = fJ-u ~ Mdi when 
fij = 0, or one in 2fij = = — /^rf, when fig = 0. Due to CP symmetry, the terms for odd = + in eq. (jS)) 
vanish. In the remaining terms, the two expansions in /is/3 and 2/1/ differ by a sum of susceptibilities with odd 
and Ud- The difference has been taken to quantify the seriousness of the Fermion sign problem [10ii| . We will use a 
somewhat different method here, as we explain below. 

The non-linear susceptibilities defined above can be written down in terms of the derivatives of Z. From the 
expression above it is clear that the derivatives with respect to the /i/ land entirely on the determinants. Now, since 
DetM = expTr logM, the first derivative gives (DetM)' = Tr (M-^M')DetM = 0iDet Af. Higher derivatives can 
be found systematically using the additional relation MM~^ = 1, which yields (M~^)' = —M~^M'M~^. Clearly, 
therefore, the derivatives of Z can be written in terms of expectation values of certain operators involving traces of 
inverses and derivatives of the Dirac operator (see the appendix for details). 

A Taylor expansion of the pressure immediately yields one for the QNS. For comparison of the notation of this 
paper to our earlier works, note that the flavour susceptibilities Xu = Xd used earlier correspond to X20 in the present 
notation; Xud, to xii! the isovector susceptibility, X3: to X20 — Xn! and the baryon susceptibility, xs = 2xo/9, to 
2(X20 + Xii)/9- The radius of convergence of the series for ^mb gives the location of the nearest critical point. The 
Taylor coefficients of X20 up to 6th order in /is/S are 

Xs=X20, xl = ^ [X40 + 2x31 + X22] , 

Xb = ^ [X60 + 4x51 + 7x42 + 4x33] , " ^ t^^" ^ ^^"^^ ^ '^^Xfi2 + 26x53 + 15x44] • (4) 

The Taylor coefficients for the off-diagonal QNS, xiii up to the same order in /is/S are 

Xb=Xi1' x| = ^ [2X31 + 2x22] , 

X%^Jl [2X51 + 8X42 + 6x33] , ^B = \\ [2X71 + 12X62 + 30x53 + 20x44] • (5) 

The coefficients of the Taylor series in 2/i/ can be obtained from eqs. 10] [SJl by flipping the sign of every NLS which 
has odd Ud- The two QNS above are in principle observable, being connected to many interesting pieces of physics 
such as fluctuations and chemical composition in heavy-ion collisions. 
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In addition, the QNS can be used to quantify the magnitude of the fermion-sign problem. Recall that the sign 
problem in the measure for the partition function in cq. ^ comes entirely from the determinant. Thus, we can Taylor 
expand the determinant to get 



DetM{n) = DetM(^o) 



l + ifj.- /io)0i + ^(a* - y"o)^02 + 



(6) 



where 0i is anti-Hcrmitcan and 02 = 0'i is Hermitean. We show later that {T/V){0i) is the number density and 
{T/V){02 + 0i0i) is the susceptibility X20- If we perform the expansion around /io — 0, then by exponentiating the 
bracket, we can see that 02 determines the width of the real part of the measure. Also, (0i0i) — {V/T)xii is the 
width of the imaginary part of the measure, up to a sign. Thus, Xii/X20 gives the ratio of the widths of the measure 
in the imaginary and real directions, and hence quantifies the magnitude of the fermion-sign problem. With a little 
care, the same argument can be extended to any fiQ. However, in that case care has to be taken to subtract the 
number density. 

We use the NLS to analyze the radius of convergence of the Taylor expansion of xb- Analysis of series expansions 
in order to extract critical behaviour is a method of long standing in statistical mechanics and lattice gauge theory 
[TgL IT^. The caveats about such analysis have clear physical meaning. First, the fact that in any infinite series, a 
few terms can be changed without changing the radius of convergence implies that one must check the series for signs 
of irregularity. In many cases these can be attributed to special features of the model (see and can be seen as 
interference from "critical" points off the real axis. We make such tests and find that the radius of convergence can 
indeed be attributed to a point on the real axis at which a divergence begins to build up. 

A. Volume dependence 

Since P, T and {/i/} are intensive variables, it is clear that each Taylor coefficient should be independent of volume 
in the thermodynamic limit. In a lattice computation, this is a non-trivial check since an arbitrary sub-expression for 
each Taylor coeflticient could diverge with volume. This divergence is canceled between terms, and with increasing 
volume, calls for ever more careful study of possible numerical inaccuracies and their elimination. A method for doing 
this was set out in 5.]. 

In order to do this systematically, we need to construct all possible sub-expressions which are independent of volume. 
At the second order, there are two sub-expressions, (02) and (0ii), which can be expressed as linear combinations of 
the two observables X20 and xii- It is clear therefore that the two operator expectation values are the basic volume 
independent expressions. The individual traces which constitute them can, and do, diverge, and it is a delicate, but 
controllable, numerical task to get a sensible answer with increasing volume. 

In we had suggested that this procedure could be generalized by notionally taking a version of QCD with larger 
number of degenerate flavours. This allows us to define a larger number of observable NLS, the maximum number 
possible at a given order iV is for Nf > N. This shows that the connected part of every fermion line disconnected 
operator must be volume independent, in agreement with other proofs, for example, in perturbation theory. 

We outline the argument at the fourth order. With four degenerate flavours at the fourth order we get — 
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The derivatives can be written easily in terms of the traces using the methods explained in the appendix- 



^4000 = Z( 01111 + 60112 + 4013 + 3022 + 04 ), ^3100 = Z{ 0iiii + 30112 



Z( 01111 -I- 20112 + 02 
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(8) 
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Finally this gives precisely the connected parts of each expectation value as the invariant - 
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(9) 



Recall that each expectation value of order two in the above expressions scales as a single power of the volume, so 
that their squares scale with two powers. This is therefore the leading divergence in the 4-th order terms, which are 
canceled to leave a result linear in volume. This first power is canceled by the factor T/V in the definition of the 
susceptibilities, giving the correct scaling of the Taylor coefficients. Such a construction generalizes. Terms of order 
2N grow as iV-th powers of the volume. The leading — 1 powers are canceled by taking connected parts, leaving 
a result linear in volume. Thus, the need for controlling numerical inaccuracies grows more severe with increasing 
order. 



B. Finite size scaling 



Consider the free energy, F{T, iJ,B\V,a,m) = F{T, fiB',V, Nt,mp,mTr), where the second expression is obtained by 
trading the bare parameters for physical quantities. Instead of the lattice spacing, a, we can specify the rho meson 
mass, nip. Then the quark mass, m, can be traded for the ratio m-^/nip since the ratio can be tuned to its observed 
value by adjusting the bare quark mass, m, all else being fixed. This process of regularizing the ultraviolet divergences 
of the theory is carried out at T = /is = 0. In other words, nip and ttItt are computed at T = /is =0. T is then 
specified by Nt and appropriate boundary conditions. 

"Small volume" finite size scaling (FSS), which we perform, establishes the minimum lattice sizes required to obtain 
the physics appropriate to the thermodynamic limit. This uses the fact that the simulations are actually carried out 
at /is = 0, below the crossover temperature for chiral symmetry breaking, Tc, by measuring the Taylor coefficients of 
a susceptibility which is expected to diverge at the critical end point. These Taylor coefficients can be written in the 
eigenbasis of the Dirac operator; for example — 

X20 ^ T 
m?p m?pV 

where |i) denotes the eigenstate of the Dirac operator with eigenvalue A^. Rewriting the coefficient, T/m^V, as 
Tm^irriTr/mp)'^ /{Vm^), one can take the factor Tm^r inside the sum, to form dimensionless combinations such as 
Tm^/X^. Then the factors m-,^/mp and Titlt^/X'^ are quark mass and lattice spacing dependent quantities, leaving 
Vrri^ to be the FSS variable at fixed T and m. Below Tc, at small values of TOttF^/'^, one is in the so-called e-region 
of physics which is dominated by the pion only. At larger values of this variable one enters the regime where the 
thermodynamic limit is approached. Experience in T = physics has shown that mT^V^^'^ «4-5 is the dividing line 
between these two regimes of physics. This is the FSS which we explore for the first time. 

Once one reaches the thermodynamic regime of sizes, further finite size scaling at the critical end point is determined 
by the divergent correlation length at the critical point. The results of ^] indicate that this is not related to the 
pion. However, as mentioned earlier, detailed "large volume" FSS of this kind, which requires resources far beyond 
those available today, will eventually furnish direct measurements of the critical indices. This we leave to the future. 
In this work we present a restricted form of this "large volume" FSS. We show that the finite size movement in the 
radius of convergence on even bigger lattices than the ones we use here is smaller than our statistical errors. 
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FIG. 2; Representing the evaluation tree for NLS. The root node is a random vector r, internal nodes stand for vectors formed 
by application of a matrix on a vector, and leaves (in black) stand for scalars. Here every vertical edge represents the application 
of on a vector (which is the costliest part of the computation), every horizontal edge stands for dotting with to give a 
scalar, and other edges for application of M' (left slanting) or M" (right slanting). The Steiner problem corresponds to finding 
the minimum cost of evaluating a given set of scalars. In the compact representation at the right, the horizontal edges have 
been contracted to bring the leaf nodes into the body of the tree and vertical edges have been contracted to a point. Note that 
the redundancy in the computation of [1 ® 2] (also written as 21) can be used as a check of the program. 

III. NUMERICAL METHODS FOR TRACES 

This section is devoted to technical matters involving the numerical evaluation of the operators 0„. We describe a 
technique which minimizes the number of matrix inversions required to determine all these operators up to a given 
maximum n. This systematic technique also allows us to identify all possible numerical tests of accuracy which come 
at negligible extra cost. We also discuss the optimization of the conjugate gradient matrix inversion. Readers with 
no interest in these details may skip this section. 

A. Evaluating traces 

The traces are evaluated through a noisy technique — 

TiA^ rt^r/ri>, (11) 

where the bar denotes averaging over random vectors r. An unbiased measurement is provided for ensembles of 
random vectors which satisfy the conditions — 



{rrYr^ - NS^pS,, (12) 

where a, (3 label the vector, i, j label the components, and may depend on the ensemble but is independent of 
either index. The matrix A consists of a product of M~^, M' and M" such that there is one for each of the 

other two. The costliest part of the evaluation is the computation of acting on a vector. While evaluating a 

certain (fixed) set of traces, we need to minimise the number of matrix inversions. 

We represent this problem through a directed graph without cycles. Each internal node of the graph contains a 
vector, and each leaf contains a scalar. The root node is the random vector r, and every other internal node is obtained 
by the action of some matrix on r. Since each vector is obtained from another by the action of either M'M^^ or 
M"M~^, the internal nodes form a binary tree. To count the cost of any operation we separate the action into a 
single followed by either M' or M" . Each internal node (except the root) gives rise to a leaf by dotting it with 
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FIG. 3: The evaluation tree for 8th order NLS requires 20 matrix inversions (numbered at the branches in the tree) and one 
extra node (in black). Since leaves are denser on the right of the tree, evaluating the inverses from right to left results in the 
use of less memory in the form of vectors to be stored. An assignment of vectors is shown (IV, 2V etc.), which uses only 4 
intermediate storage vectors. Bidirectional arrows connect quantities that serve as checks of the program. 

r^. Since this action is the stochastic evaluation of a trace, several internal nodes (differing by cyclic permutations of 
the operations) may connect to each leaf. A representation of the evaluation tree up to level 2 is given in Figure [3 A 
compressed representation, also shown in Figure [3 is obtained by collapsing together the nodes containing a vector 
V and the vector M~^v, by writing 1 for each Af, 2 for each M" and never writing M^^. 

Our problem is — given a set of target leaf nodes, find the path on the tree using the minimum number of M^^ 
which evaluates these leaves. This is one version of a problem known in the computer science literature as the group 
Steiner problem on directed graphs, for which a solution has been given recently ^19J. The computer science interest 
arises from the fact that the general problem is known to be in the class of NP complete problems. 

Since our problem size is small, we do not use the general algorithm, but a heuristic which shares the idea of 
"bunching" with that algorithm, uses the structure specific to this problem, and is easily implemented in wetware. 
Given the set of matrix elements, write down all strings which have to be generated, and group each with all its 
cyclic permutations. For each set of strings of length £, starting from the largest, enumerate the suffixes of length 
£ — 1. Choose that representation in which the suffixes are common at the nearest possible level. This is the central 
heuristic — bunch the strings by largest suffixes. Build this back all the way to the empty string, backtracking only 
if there is insufficient bunching close to the root. This gives an evaluation tree. Run over all such trees and find the 
ones with lowest cost. 

For the 4th order problem this yields: 
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FIG. 4: The evolution of the norm of the residual, |r„p against the number of conjugate gradient iterations, n. This figure 
shows the computation on a typical gauge configuration at T = 1.05Tc using one random vector for the full evaluation tree of 
the 8th order NLS. In Method 1 n increases as one proceeds down the tree. In Method 2 n remains almost constant (see the 
thick line). 

(3-1) : 111, 
(2 • 1®2) : 112, 121, 211 
(4-1) : im, 

(13) 

The underlined suffixes are used in the evaluation. The tree is identified from bottom up, but evaluated top down. A 
solution for the 8th order problem is given in Figure |21 

There seems to be a large degree of redundancy in the most efficient evaluation. This can be used to minimize 
the number of internal nodes not connected to the target leaf set. Even after this is done, there seems to be further 
degeneracy, and this can be utilised to minimize the breadth of the tree, since that determines the number of vectors 
needed in the evaluation. 

Sometimes only one of the two branches is needed in the evaluation. Since both branches can be generated at 
nearly the same cost as one, the other branch can then be used to check the numerical accuracy of the procedure. 
Such possible checks are marked by bi-directional arrows in Figure |31 

B. Inverting the fermion matrix 

The fermion matrix inversion is done through a conjugate gradient method. The most time consuming part of this 
problem is the matrix multiplication M'' Mx. A huge gain could be obtained if the multiplication could be speeded 
up in any way. Each matrix multiplication involves V diagonal terms and 2dV off-diagonal terms (where V is the 
number of sites, and d the number of dimensions). Each diagonal term involves Nc multiplications. Each off-diagonal 
term involves 2N^ operations {N^ multiplications and additions). Hence the total operation count for the double 
matrix multiplication is VNc{l + AdNc). 

One consequence of this is that 'improved' methods such as preconditioning or deflation are not of much use unless 
the number of matrix multiplications can be cut down drastically. For example, in our application where the same 
matrix is applied repeatedly on many different vectors, one might expect deflation to work wonders po|. However, 
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FIG. 5: Results of inversion with e = 10"'' with method 1 (fuU lines) and with e = 10~^ with method 2 (dotted lines). The 
latter are significantly different for {ipip) and Xir, but not for the traces &„■ Inversions with e — 10~^ with method 2 are 
indistinguishable from the full lines on the scale of this figure. The panel on the right displays the accuracy of the inversions as 
one proceeds down the evaluation tree. The significance level is defined as the difference of two evaluations of the same trace 
divided by the error in the difference. 



this requires two double matrix multiplications per step and hence actually performs worse. 

The simplest solution, and the one that seems to work best, is to tune the stopping criterion. In Figure 21 we show 
the performance of two different criteria — 

1. Method 1 is to require the norm of the residual, |r„p < eV . Since each inversion increases the norm of the 
right-hand side, through the 20 inversions required for evaluating all the traces, the CG has to take increasingly 
larger number of steps, n, for convergence. 

2. Method 2 is to require the norm of the residual, |r„p < e|rop. As shown in the figure, this seems to keep n 
almost fixed as one proceeds through the evaluation tree. 

There is really no a priori reason for the two stopping criteria to give equally accurate results with the same value 
of e. However, we see that for e < lO"** the results do not differ for the 0„'s (they do for (ip^p)), while method 2 
takes less CPU time. To have some understanding of this, we note that since any matrix is diagonalised by a bi- 
unitary transformation, we can write M = U^DV, where D is diagonal, and U and V are unitary matrices. Similarly, 
M" = WWX, giving 

d" 

Ti- M-Hl" = Ti D-^UW^D"XV^ ^ ^ -f-{UW^)^j{XV^)j^ (14) 
The triangle inequality can then be used to bound 



\ttm-^m"\ < 22 
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Since every element of an unitary matrix is bounded by unity, we get the second bound above by inserting unity for 
these matrix elements. However, this bound is too rough. If it is (nearly) saturated then the trace would depend 
on the stopping criterion as strongly as {4>^) — TrM^^, contrary to observation. Therefore it seems likely that the 
structure of the unitary matrices plays a role in reducing the bound. In particular, we note that the low eigenvalues 
of M, da < p can be killed only by eigenvalues of M". This can happen if the unitary matrices WU"^ and XV'' are 
nearly block-diagonal in this subspace (they are exactly block-diagonal in the vacuum or in the presence of stationary 
gauge fields, when the Dirac operator becomes separable). In other words, the observed lack of sensitivity to small 
eigenvalues of M can be understood if every near-zero mode of M is also a near-zero mode of M" . 
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FIG. 6: The Wilson line susceptibility, xl, shown on the left as a function of the bare coupling, l3 on lattice sizes of 8"^ 
(triangles), 12"^ (pentagons), 16^ (boxes) and 24^ (circles). It is clear that /3 — 5.2875 does not locate the exact position of 
the peak on the larger volumes. This is also indicated in the figure on the right, which depicts X-C a-t a fixed l3 = 5.2875 as a 
function of V. 



IV. LATTICE SIMULATIONS AND RESULTS 



A. Parameters and simulations 
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TABLE I: The couplings (3 and bare quark masses m at which the simulations were performed. The temperature scale, T/Tc, 
is given in the MS scheme, with error assignment as discussed in the text. The bare quark mass at all temperatures correspond 
to m = O.lTc. The statistics collected in each simulation, Nstat, is given in terms of the longest autocorrelation time, rmax, 
measured in that simulation in MD time units. The trajectory lengths are 1 MD time unit on the 4x8'^ lattice and scaled up 
in proportion to the length of the spatial lattice, i.e., growing to 3 MD time units on the 4 x 24'' lattice. 



All the simulations reported in this paper were carried out at the lattice spacing a = 1/4T, corresponding to lattices 
with temporal extent Nt = 4. We used the R-algorithm in the simulations and chose the conjugate gradient stopping 
criterion to be Method 1 of Section UTTl with e = 10^^. The finite temperature crossover at this cutoff has been 
determined to be at the coupling /? = 5.2875 ± 0.0025 |23|, which therefore corresponds to the cross over temperature 
Tc at zero /x. We have improved the bound on the cross-over coupling, as we report later in this section. 

We have set the scale by using plaquette values to define a renormalized coupling |22l | in a series of exploratory 
computations on small lattices (16'') at T = 0, using trajectory lengths of one molecular dynamics (MD) unit. The 
typical statistics utilised for scale setting was about 500 such trajectories, of which the initial 200 were discarded for 
thermalization. During the scale setting the input bare quark mass was varied until a self-consistent determination 
of the scale showed that m/T^ = 0.1 within statistical errors. A two-loop analysis of the scale shows several sources 
of uncertainty. Statistical uncertainties can be easily made small, and with our statistics they are of the order of 3~6 
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FIG. 7: The panel on the left shows operator expectation values (0n) as functions of the volume at T = O.TSTc. The panel on 
the right shows that the connected parts {022)0, and (0222)0 at 0.75Tc scale with the correct power of V . 
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FIG. 8: A comparison of Gaussian noise (boxes) and Z2 noise (circles) for the evaluation of the NLS X40 (left panel) and xeo 
(right panel) at 0.75Tc on a 4 x 16^ lattice. Gaussian noise requires less vectors for stable evaluation of these quantities. 



parts per thousand. Harder to control are the differences between various renormalization schemes and the scatter 
due to errors in the estimate of the critical coupling. The latter induces an error of slightly over 1% in the scales. 
Differences between different renormalization schemes indicate that a = l/iTc is too coarse for 2-loop formula; to 
work with high precision. This problem is also known in the quenched theory, where it was seen that a higher order 
correction was needed to achieve scaling 22,3. Here we quantify this uncertainty by taking the scatter between the 
scale determination through the E, V and MS schemes. The uncertainty defined in this manner is within 2% of 
the value in the MS scheme for T within 25% of Tc, but grows rapidly beyond, upto about 10% at the value of (3c 
determined a,t Nt — 8 for m/Tc = 0.1. Since the critical end point lies close to Tc, we decided not to do an extensive 
set of T = simulations to control the scale further by extracting higher order terms in the /3-function. The combined 
statistical and systematic errors in the determination of T/T^ are shown in Table 

Hadron masses have been determined earlier at the critical coupHng It was found that m,r/mp — 0.31 ± 0.01, 
which is somewhat larger than in nature. Tc could be given in physical units by noting that rrip/Tc = 5.4 ± 0.2. 
Further evidence that lattice spacings are still fairly coarse can be seen by noting that the nucleon is too heavy, since 
rriN/iTip = 1.8 ± 0.02. Pushing towards the continuum will be a necessary future step. 

The simulation parameters and statistics are shown in Table The MD trajectory was chopped into time steps of 
length 0.01, and the full trajectory was chosen to be 1 time unit for the 4x8'^ runs. As the spatial size, Ns of the 
lattice changed, the trajectory length was changed in proportion to N^. An earlier study had shown such a tuning 
of the trajectory length to be useful in controlling the growth of autocorrelations We found that the number of 
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FIG. 9: Estimates of Xii/^^ as functions of the number of vectors for O.TSTc on a 4 x 24'^ lattice. Shown are the estimates 
using Gaussian noise and CG stopping criterion of 10~* (circles) and 10~^ (boxes) in Method I and Zi noise with CG stopping 
criterion of lO"** in Method I (triangles). are chosen to be multiples of 10, but some of the points are displaced slightly in 
the x-direction for improved visibility. The panel on the right shows x^o as a function of Nv for n = 4 (circles) , 6 (boxes) and 
8 (triangles). This data was collected at O.TSTc on a 4 x 10"^ lattice with Gaussian noise and eca ~ 10~*. 



iterations in the conjugate gradient algorithm remained unchanged as the spatial volume increased at fixed T. 

The pseudo-critical point is usually identified by the Wilson line susceptibility, yt , the chiral susceptibility, Xm, 
which is the second derivative of the free energy with respect to the quark mass |25| , or autocorrelations of various 
measurements. All of these are expected to peak at T^. If the point is critical, then one expects these quantities to 
grow as a power of the volume, V . If, on the other hand, there is a cross over, then these peak heights should saturate 
at large enough V . Also, different quantities could peak at somewhat different temperatures even for F — > oo (Tsj . 

In FigureElwe show a measurement of xl at the couplings we have used. The drop in the peak of xl on the largest 
V is clear indication that /? = 5.2875 is not the location of the peak in the thermodynamic limit. In order to bracket 
the position of the peak we have run separate simulations on 4 x 12^^ and 4 x IG'^ lattices at nearby couplings with 
bare quark mass fixed at ma = 0.025. Collecting statistics from about 30 to 50 independent configurations at each 
coupling, we find that the finite volume shift in the peak is less than A/3 = 0.00125 for Ng = 16. This corresponds to 
an error of 1% or less on T^- 

Such a shift is not of significant concern for our study since the dominant systematic error in the scale setting comes 
from the fact that the lattice spacings for A't = 4 in the temperature range of interest are too coarse for 2-loop scaling 
to work — improving the estimate of beyond the present 1% level of scale does not improve the total systematic 
error by 1%. Of course, as a consequence, these data cannot yet be used to check whether there is a phase transition 
or a cross over in the /i = theory. This is an interesting but separate problem "26", '27] . It might be useful in future 
to explore a narrow region of fi near the peak in greater detail to track the growth of xl with V and give a definitive 
answer to this question. 

Several local gauge operators as well as quark operators such as the chiral condensate were measured once per 
trajectory. The longest autocorrelation, Tmax, in these measurements was used to set the scale of autocorrelations. It 
was seen that Tmax, when measured in MD time units, was roughly independent of Ns at fixed T, as expected from 
the scaling of the trajectory length (2^. As a consequence of these scalings, the time taken for generating one new 
configuration scaled approximately as iV^ at fixed T in thermal equilibrium. 



B. Noise and measurement 



Full control over errors in determining the non-linear susceptibilities is crucial for the job of extrapolation to non- 
zero \i. The major source of numerical errors is the numerical cancellation of divergences which go as powers of the 
volume in order to get a finite result. One part of this program is the identification of primitive volume independent 
quantities — the connected parts of fermion-loop disconnected operators, some of which were enumerated in Section 
III Al In this section we complete the program by showing that some choice of the noise ensemble works better than 
others p3 |. and that the number of vectors is the crucial parameters, not the choice of the conjugate gradient stopping 
criterion. 
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FIG. 11: A fourth order QNS and some connected parts on lattice sizes 4 x 12"^ (triangles), 4 x 16"' (circles) and 4 x 24^ (boxes). 
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FIG. 12: Some connected parts on lattice sizes 4 x 16^ (circles) and 4 x 24^ (boxes). The first panel shows (T/V^)(022). The 
second panel shows a model cross-over (full line) with first (dashed) and second (dotted) derivatives. The third panel shows 
(T/V){0222)c. Note the similarities between data and model. 



In Figure[7|we show the easiest part of this program — the scaling with volume of the expectation values [T /V){On) ■ 
These are fairly easy to handle because the operators have no disconnected pieces and there arc no divergences to 
cancel. We find that the values and the error estimates are very stable against the choice of the conjugate gradient 
stopping criterion, in method I we can change e from 10~^ to 10^^ without making a change to these quantities. 
There may be a little volume dependence in going from TV^ = 8 to A^s = 12 lattices, but this saturates for Ng between 
12 and 24. The errors become larger as n increases. 

We found that, in agreement with expectations, the disconnected pieces of traces such as {Qij...) grow with volume 
as a power equal to the number of fermion-line disconnected pieces in the operator. However, as shown in Figure 
[3 the connected part grows roughly linearly with volume. The residual volume dependence is small, and the fact 
that {T /V){&222)c decreases with V shows that the cancellation of the leading divergent pieces has been performed 
satisfactorily. These results have been obtained with 100 noise vectors drawn from a Gaussian ensemble. 

We have made further tests of the effect of the choice of the noise ensemble, as shown in Figure |H1 It seems that 
it is preferable to use Gaussian noise rather than Z2 noise, since the former is more stable against changes in the 
number of random vectors used. Since quite the opposite conclusions have been presented in the literature, albeit 
for different measurements [28[, we investigated this a little further. It seems that indeed Z2 noise performs better 
for single traces, as in the measurements of (0„). However, the measurement of products of traces such as (0^...) is 
cleaner with Gaussian noise. 

The two final issues to be resolved are the role of the accuracy in the conjugate gradient, and the number of 
vectors required for measurements. In Figure |51 we show that the choice of the CG stopping criterion, ecGj does not 
significantly affect the choice of Ny. It is also clear from this figure that larger values of are required for the 
higher order measurements. In particular, Ny = 100 seems to suffice for measurements up to the 4-th order, whereas 
Ny = 200 is required beyond that for measurements up to the 8-th order, in the hadronic phase. In fact, iV„ also 
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FIG. 13: The radius of convergence of the series expansion for X20- The first panel shows rg at T/Tc — 0.9 (pentagons), 0.95 
(circles) and 1 (triangles). The second panel shows results at T/Tc = 0.95 for (circles) and (boxes). 




FIG. 14: We show the dimensionless series coefficients, Xs^^""^ (^q. BJj obtained on 4 x 24^ lattices, for n = (pentagons), 
n = 2 (circles) and n = 4 (boxes). The n — 2 coefficients have been divided by 2 and the n = 4, by 12, in order to bring them 
all to roughly the same scale. 



depends on the lattice size. On smaller lattices N.^, = 200 seem to suffices, but this grows with lattice size and on the 
4 X 24"' lattices one needs Ny — 500 below Tc- Above Tc, Ny ~ 100 seems to suffice for all the measurements and 
lattice sizes. 
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FIG. 15: X2o(Ms)/r^ obtained using merely two terms in the Taylor expansion already shows the qualitative change above 
T = 0.9Tc. The 66% confidence intervals for the extrapolation are shown for Tc (full lines), O.OSTc (dash-dotted) and O.QTc 
(dashed) from computations on a 4 x 24'^ lattice. 




FIG. 16: We show the radii of convergence as a function of the order of the expansion at T = 0.95Tc on a 4 x 8'^ lattice (circles) 
and a 4 X 24^ lattice (boxes). The panel on the left shows pn and the one on the right is for Tn- 



C. Susceptibilities and radii of convergence 



In this subsection we present results on the temperature dependence of the linear and non-linear quark number 
susceptibilities. Comparable data on the linear QNS in the high temperature phase were presented in [29j . In figure lTUI 
we present data for Xi^jT"^ and Xn/T^. Some volume dependence is visible in the immediate vicinity of Tj,. The high 
temperature behaviour is compatible with the predictions of [23, 1^ • The results are also completely compatible with 
earlier results in [23| after correcting for a division by an extra factor of (T/F) for Xii/^^^ reported there. Comparison 
with the recent results of |32| are harder to perform since the actions are different. As a result, such a comparison 
can be meaningfully performed only after taking the continuum limit. 
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There is a clear qualitative difference in xii above and below Tc. Above Tc it is significantly easier to measure, is 
small and compatible with zero. Below, and in the vicinity of T^, it is possibly larger, but much harder to measure 
because of fluctuations. The ratio Xii/X2, taken at fis — 0, quantifies the magnitude of the fermion sign problem 
at small /is , as discussed in Section ^1 This ratio is also shown in Figure ^| It is clear that the fermion sign 
problem is negligible at higher temperatures. The success of modern extrapolation methods at finite /i ^ .lOj is, 
at least partially, due to this. Due to the rapid increase in phase fluctuations on lowering T, visible in FigureEB any 
extrapolation technique would be hard to use significantly below the lowest temperatures we have used. Since xii 
depends strongly on the pion mass, we expect phase fluctuations to pose greater difficulties as one approaches the 
physical pion mass. Conversely, the sign problem would be easier to deal with in simulations at higher pion masses, 
such as • 

We show the fourth order NLS, X40: in Figure^Jfor our three largest lattice volumes. Away from the critical region 
the volume dependence is seen to be negligible. In the critical region, however, Xio exhibits behaviour similar to that 
of XL (Figure 1^1). The susceptibility peaks near Tc and shows pseudo-critical behaviour, since the near-critical region 
visibly shrinks with increasing volume. Since X40 decreases with V at fixed /? — 5.2875, we again have evidence that 
the cross-over coupling is near this, but not exactly here. 

Further insight into the nature of this peak comes from examining the other NLS at this order. It turns out that 
X22 also has a similar sharp peak, whereas X31 is much smoother. These lead us to examine the connected parts 
which contribute at this order (see eq. EJ. We found that {T/V){04), shown in Figure ITTI seems to behave almost 
as a (pseudo) order parameter for the crossover at T^. In this it behaves like X2o/T'^ and, indeed, like all the single 
traces {T/V){0n) that we have examined. The peaking behaviour is essentially due to a similar peak in {T/V){022)c, 
which we have displayed in Figure [T^ 

Note that 02 is a composite bosonic operator. It may be possible to write down effective long-distance theories in 
which it is treated as a field operator whose expectation value shows the correct cross over behaviour. In that case 
(022)c would be the susceptibility of this field, and being proportional to the temperature derivative of (02), would 
peak, as observed. In fact, one can push this idea further, and examine the expectation values (0222)c. In such an 
effective theory one would find this to be proportional to the next derivative of (02). Then the T-dependence of these 
quantities at fif — would have the shapes shown in the second panel of Figure [T^ In fact, the data, shown in the 
final panel of Figure [T^ does show the same qualitative behaviour. Thus, the cross-over at Tc can be probed by any 
of these connected parts. 

There is significant volume dependence in the close vicinity of Tc, although for T < 0.95Tc one is clearly outside 
the pseudo-critical region, and the volume dependence is weak. Within the pseudo-critical region the peak of (022) c 
is seen to have the same qualitative behaviour as the peak of xl shown in Figure including similar evidence for 
finite volume shift in the cross over coupling. There are consequent effects on the higher order connected parts, as 
shown in Figure [T^ 

With the NLS at hand, we construct the Taylor expansion of the diagonal QNS, X2o(M-b): where the coefficients 
are given in eq. Radu of convergence can be obtained for the series expansion f{x) — J2 hriX^^ using two 

definitions — 



Pn = 



/2 



l/2ri 



and 



/2n 



2n+2 



1/2 

(16) 



At large n both definitions should be equivalent, although for many series it is known that the latter gives somewhat 
smoother approach to the limit of n — > 00. Also, at a critical point, one may extrapolate the latter to the limit n — > 00 
by a fit to a 1/n behaviour 16]. We apply these definitions to the series for the baryon number susceptibility (see 
above eq. . 

In Figure [T^ we display the radii, r4 and rg, at several different temperatures as a function of N^: the volume 
dependence is strong. For 0.95 < T/Tc < 1, the radii roughly agree with the estimate of the critical end point in 
for Ns comparable to those used in that study. However, as the volume increases, the radius of convergence drops. 
The magnitude of the finite volume effects is parametrized by the pion Compton wavelength, and for m^A^^ > 5-6 
one finds that these effects saturate. This is in accord with the discussion in Section Hi Bl Similar saturation of finite 
size effects has been seen in other measurements in the chiral sector at similar values of nij^Ns |T^ . 

At temperatures below 0.95rc the radius of convergence falls quite dramatically. However, examination of the series 
coefficients shows that in this lower range of temperatures the 6th order coefficient, x%7 is negative (see. Figure 
Thus the radius of convergence does not indicate a divergence. Rather, it shows that at this temperature X2o(M-b) < 0, 
thus violating thermodynamic stability, and therefore indicating that higher order terms in the expansion become 
necessary. This is consistent with the increase in the ratio X11/X20 at these T (see Figure At temperatures 

between Tc and 0.95Tc, the character of the expansion is quite different, with all the Taylor coefficients being positive. 
As discussed earlier, when this is true the series diverges on the real fiB axis. The qualitative difference in these 
two expansions is already visible at low order, and is illustrated in Figure [T^ Note that a figure like this illustrates 
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TABLE 11: Summary of critical end point estimates — the lattice spacing is a = l/4r. Ns is the spatial size of the lattice 
and NsiriTr is the size in units of the pion Compton wavelength, evaluated for T = n = 0. The ratio mTr/niK sets the scale of 
the strange quark mass. As the mass scales indicate, the lattice spacings and u and d quark masses of and this work are 
comparable. The numbers in brackets indicate the errors on the least significant digit. Since no estimates of the critical end 
point are quoted in l9lll0|l. the comparison is more complicated, and is described in the text. 

whether or not there is a divergence, but since it shows the sum of a finite number of terms, it cannot show anything 
special happenning at the radius of convergence. 

In Figure^lwe show the variation of the two different definitions of the radius of convergence, pn and r„, with the 
order n. On our smallest lattice, 4x8^, the successive estimates diverge, and the two sets of estimates are compatible 
at the 2-a level. After the cross over to large volume behaviour, the estimates stabilize with order, being nearly 
independent of n. Also p„ and r„ are in close agreement with each other. The radius of convergence extrapolated to 
ah orders at = 0.95Tc gives /x^/T^ = 1.1 ± 0.2. 

If we assume that the end point is in the Ising universality class, then the finite volume shift in the end point is 
given by 

,^{V) = ,-{^) + ^, ^.e. ^.-iV)-p-{V') = ^-^, (17) 

where ^ = 5 is the 3-d Ising magnetic exponent and V = N^. Using the data exhibited in Figure [TSl for Ng — 16 
and 24, which are both above the crossover, it is easy to check that the finite volume shift in fi^ is bounded by the 
statistical error in the estimate for Ng = 24. In view of this, we quote the estimate and its error obtained for the 
Ns = 24 lattice as the estimate in the thermodynamic limit. 

V. SUMMARY AND DISCUSSION 

The main result we present in this paper is the strong volume dependence of location of the critical region, leading to 
a substantially smaller estimate of /Xg on large volumes than before. We began by constructing the Taylor expansion 
for the quark number susceptibility, X2o{T, pb), for Nf = 2 QCD at for several large volumes. The leading (/is 
independent ) te rm of the series is the quark number susceptibility (QNS) which has received extensive attention 
recently p9ll30l [Sllls^ IsM Is ^ . The remaining Taylor coefficients involve the non- linear susceptibilities (NLS) which 
were defined in !5|. We have taken the series up to the 6th order term, which involves 8th order NLS. A systematic 
and efficient procedure for generating and computing the quark number susceptibilities at any order was presented in 
Sections im and ITTIl 

The main thrust of this work is to approach the large volume limit. From the point of view of thermodynamics, 
this limit has to be taken in order to decide whether rapid changes in certain quantities and peaks in others are due to 
critical behaviour, cross over or a first order phase transition. For the fj.B — cross over at Tc, good progress has been 
made |i26> |. However, due to questions such as the magnitude of the shift in the critical coupling, and the absence of 
evidence for peaking of different susceptibilities at slightly different points, an unqualified answer is still not available. 
Note, however, the remarkably high accuracy in the measurement of Tc that is possible even without answering this 
question. The situation is similar at critical points. It is usually much easier to identify the critical point (assuming 
it is critical) than it is to "prove" criticality through extensive analysis. In this work we have shown that there is a 
minimum lattice size that one should use in order to get a reliable estimate of the position of the critical end point. 
Further work involved in "proving" criticality by studying detailed finite size effects, with substantially larger lattice 
volumes and extracting critical indices, is beyond our present computational capability and is left for the future. 

Taking the large volume limit is also necessary to resolve the small eigenvalues of the Dirac operator, and thus 
control the chiral behaviour required to see the critical end point. Since we use large lattices as well as small quark 
masses, there are new technical problems to be solved (see Section FlV B|l . Using the methods developed here, the 
time required to compute the expansion to 8th order is between 1 to 3 times that required to generate an independent 
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FIG. 17: The phase diagram and the freeze-out curve superimposed. The filled circle denotes the estimate of the critical point 
which was obtained in this work. The open circle is an earlier estimate from using smaller lattices and nearly the same 
quark mass as in this work (see the text for a discussion of their comparability). The freezeout curve has been converted 
for this figure using a value of Tc appropriate to our computation. A scale of the CM energy per nucleon, \fS, has been marked 
on the freezeout curve. 

configuration well away from the finite temperature cross over at Tc. Closer to Tc the time required for computing 
the Taylor expansion drops rapidly to 1/10 of the time required to generate a statistically independent configuration. 
Since one never actually computes anything at the critical point, there is no critical slowing down. 

Away from Tc we found little volume dependence for the QNS and the NLS. This extends our earlier results, which 
were confined to T > Tc |29l |. The off-diagonal quark number susceptibilities, xiij reasonably large (but noisy) 
below Tc and become small in the high temperature phase, where they are in rough agreement with the results of 
a perturbative computation in (3Q | . The ratio of these two susceptibilities is a direct estimate of the severity of the 
sign problem (see eq. 01. We found that the sign problem is very severe at small T, and acute even at T « 0.85Tc 
(see Figure [TU)) . This manifests itself in the Taylor expansion by requiring that many terms in the series be summed 
in order to get physically sensible results (see Figure [T51 and the discussion in Section llV C|l . This aspect of the sign 
problem has been investigated in [s^ and named the "silver blaze problem" . 

Near Tc, on the other hand, volume dependences were significant. This is, of course, necessary for a system to 
become critical. However, we presented additional evidence that there is no criticality at /is = b y ex amining the 
volume dependence of the maximum of xl (Figure IHJi. This adds to the growing body of evidence [ll,|2^|23| that 
the critical point moves out to finite /is when the pion mass (or equivalently, the quark mass) is finite. Finally, 
the analysis of the series expansion of X20 /T"^ showed clear evidence of a breakdown and consequent criticality for 
T^ K, 0.95 and fj,g/T^ « 1.1 on large lattices (see Table Hi) for the error bars). We saw a crossover from small to 
large volume behaviour for lattice sizes A^gW^ « 6 (see Figure [T5|) . For the pion mass used here, that corresponds to 
A^s ~ 14. Similar finite size effects have been reported in 18!], and ascribed to the lack of energy resolution for the 
small eigenvalues of the Dirac operator when working at small volumes. 

In Table Hll we have collected all estimates of the critical end point in recent lattice computations, along with other 
relevant scales. It is interesting that the nucleon mass is still very large compared with the real world. While this mass 
would change on tuning the quark mass, the ratio mi^/mp will not decrease significantly towards the continuum value 
except by changing the lattice spacing. Since the nucleon mass is important for determining the value of the QNS X20 
below Tc, it is clear that going towards smaller lattice spacing is equally important for getting a reliable estimate of 
the critical end point of QCD. Unfortunately for computations with improved actions, there is as yet no evidence that 
mAr/wp takes on its physical value already at coarse lattice spacing. We emphasize that our method is capable of 
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handling the larger lattices required to take the continuum limit; the only constraint is the usual polynomial growth 
in the required CPU time as a function of N^. . From a comparison in Table [H] of our results with those obtained in 
0,01, it appears that the strange quark in seems to play only a small role in the determination of the scales. 

Instead of an estimate of the critical end point, ^ quotes an estimate of the pseudo-critical line (3{ii), in the 
expectation that the critical end point lies on this. The end point estimate of 6] is indeed comparable. On comparable 
lattice volumes we found that our method of estimating the critical end point gives results which are comparable to 
0,1^. Our true estimate of the end point, quoted in Table UTI however comes from computations on larger volumes, 
as explained above. 

The estimate of the end point in is indeed carried out on large volumes in units of the pion Compton wavelength. 
However the pion mass used there is twice as large as here, and hence far from the physical value. The compatibility 
of the results of P| and ^3 can then be ascribed to a fortuitous cancellation of volume dependence and the pion 
mass effects. Knowing that in these two computations differed by a factor of two, this conclusion could have been 
reached earlier by noting that the value of /i^ /T^ drops by a factor of two when is decreased by a further factor 
of two [12. 

This is an appropriate place to list caveats. By identifying the radius of convergence with the location of the critical 
point, we have found a fairly broad region over which indications of criticality can be seen. Strictly, the radius of 
convergence is just the lower bound for the location, and higher order terms should be examined in future. The 
extent of this critical region in the /is direction is marked out by the filled circle in Figure [T7I The fact that present 
computations, ours included, are performed on coarse lattices with lattice spacing a — 1/4T has several implications. 
The most important is the question of scale. At such lattice spacings the ratio mjss jrfip is too large. This implies that 
the value of changes by a large factor if one sets the scale by the nucleon mass rather than by the rho mass. In 
effect this means that one has to go to finer lattice spacings to get the scale. We have set the scale by the somewhat 
more stable method of using t\-^ identified from a running coupling on the lattice [23. Is^ . However, one needs to 
go to smaller lattice spacing in order to stabilize this estimate (see Section llV Al for details). It has also been pointed 
out in [3 that the continuum limit has to be taken in order to remove an ambiguity in the prescription for putting 
chemical potential on the lattice. 

We further note that direct proof of criticality through the extraction of the critical indices at the end point would 
require significantly longer series expansions than we have used here. It would be useful to convert the algorithms for 
generating the coefficients into "compilers" so that this job can be done automatically. In the interesting region near 
Tc the measurements take a small fraction of the time needed to generate the configurations. As a result it seems 
plausible that such extensions of the present work can be undertaken in the near future. 

Nevertheless, it is interesting to speculate on questions whose answers are beyond our ability to compute at present. 
Since /z^ /T^ is expected to decrease with ttItt , our estimate of the end point puts an upper bound to this quantity. 
Finite lattice spacing errors were estimated in [^ and can be controlled in computations in the near future. This 
upper bound on the estimate of the critical end point we obtain suggests that an experimental search for the end 
point is feasible today. 

We have illustrated this in Figure [T7I bv plotting our end point estimate along with a freezeout curve determined 
in |3It| . This freezeout curve has been converted to T/Tc units using a value of appropriate to our computation. 
Since this curve is obtained by treating a resonance gas as ideal, we expect this freezeout curve to lie below the line 
of transitions and the critical point, as it indeed does. We have marked out on the freezeout curve values of the CM 
energy, v^, per nucleon needed to reach that point on the curve. The fact that the heavy-ion experiments at CERN 
with = 18 GeV did not see any sign of a nearby end point is a confirmation that the pion mass in the real world 
is lower than that used here. An energy scan at the RHIC should be able to locate the critical end point of QCD. 

Acknowledgements: We would like to thank Simon Hands and Owe Phillipsen for discussions during the program 
"QCD at finite density: from lattices to stars" of INT Seattle. SG would also like to acknowledge Jaikumar Rad- 
hakrishnan for developing the connection between the problem of efficient trace evaluation and the Steiner problem, 
Keh-Fei Liu for a discussion on the properties of Zi noise, and Edwin Laermann for discussions during a visit to 
the University of Bielefeld funded by the French-German Graduate school funded by DFG under grant No. GRK 
881/1-04. RVG would like to thank the Alexander von Humboldt Foundation for financial support and the members 
of the Theoretical Physics Department of Bielefeld university for their kind hospitality. This computation was carried 
out on the Indian Lattice Gauge Theory Initiative's CRAY XI at the Tata Institute of Fundamental Research. It is 
a pleasure to thank Ajay Salve for his administrative support on the Cray. 



21 



APPENDIX A: THE SUSCEPTIBILITIES 
1. The susceptibilities 



There are two steps to writing the NLS. First the derivatives of the free energy (or log.^) are expressed in terms of 
the derivatives of Z. Second, the derivatives of Z are expressed in terms of products of traces (the quark operators). 
The first step is to take the derivatives — 



a"P(r,{M}) 

9/Uj9/Uj • • • 



T\ S"«'+"<'logZ(T,{/i}) 



(Al) 



This is easily accomplished in any computer algebra system; for example, in Matliematic;a using the program 
fragment — chi [nu_,nd._] := D [Log [Z [muu.mud] ] , {muu.nu} , {mud.nd}] augmented by substitution rules for set- 
ting Xn„,nd = foi' odd Uu + Ud, and implementing the symmetry n„ <-> Ud- The results are written out below. Since 

Z is a moment generating function, log ^ is a cumulant generating function. Statistics mavens will therefore recognize 
the expressions below as the definitions of cumulants in terms of moments. At the leading order we have 



Xio 



T\ Z^ 
V Z ' 



(A2) 



Since this is zero, we shall use the relation Zio = to simplify successive derivatives. At the second order we find 

T\ Z20 fT\ Z\\ 



X20 



Xii = 



V z 



(A3) 



At the third order we have we obtain the relations Z30 = Z^x = which can be used to simplify the later derivatives. 
At fourth order we have 
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X31 
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At the fifth order CP symmetry allows us to write Z50 = Zn = Z32 = 0. The 6th order susceptibilities are — 
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The 7th order susceptibilities give the further relations Z-^q = Zgi = Z52 = Z43 = 0. At the 8th order we obtain the 
susceptibilities 
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The second step is to write the derivatives of Z in terms of products of traces. A diagrammatic method for this 
has been given before [s^. We put down a second method suited for symbohc manipulations, which, however, is 
recursive. The two basic identities are Z\q = Z(0i) and 0^ = 0n+i- These give first the low order derivatives — 

(A7) 



Z2o = Z{0n + 02j, Zn^Z{&n^ 
Here the notation stands for the product, 0i0j • • • 0;. At the fourth order we get- 

Z40 = ^( 01111+ 60112 +4013 + 3022 + 04 



^31 - ^^01111 + 30112 + 013), 
Z22 = ^{^01111 +20112 + 022 

At the 6th order the derivatives are 
Zgo 



(A8) 



^(^0111111 + 15011112 + 2001113 + 4501122 + 150114 + 600123 
+6015 + 150222 + 15024 + 10033 + 06 
Z51 = Z<^0iiiiii + 10011112 + 1001113 + 1501122 +50114 + 100123 
Z42 = Z<^0iiiiii +7011112 +401113 +901122 + 0114 +40123 + 30222 + 024 

^33 = ^<^0111111 +6011112 + 201113 + 901122 + 60123 + 033^ (A9) 

Finally we write down some of the eighth order derivatives — 
Zso = Z( 011111111 + 2801111112 + 560111113 + 70011114 + 21O0iiii22 + 560iii5 + 56O0iii23 + 280ii6 
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+420011222 + 42001124 + 28001133 + 8017 + 84001223 + 1680125 + 2800134 + 10502222 
+2100224 + 2800233 + 28026 + 56035 + 35044 + 08 



^71 = Z{ 011111111 + 2101111112 + 350111113 + 1050111122 + 35011114 + 21O011123 + 2101115 + 1O5011222 

+7001133 + 10501124 + 70116 + 10501223 + 350134 + 210125 + 017^ 
^62 = Z{ 011111111 + 1601111112 + 200111113 + 600111122 + 15011114 + 80011123 + 601115 + 60011222 + 1001133 

+3001124 + 0116 + 60125 + 6001223 + 1502222 + 150224 + 100233 + 026 
Z53 = 2'( 011111111 + 1301111112 + 110111113 + 450111122 + 5011114 + 50011123 + 01115 + 45011222 + 1001133 
fl50ii24 + 450 1223 + 30125 + 50134 + 100233 + 035 
Z{ 011111111 + 1201111112 + 80111113 + 420111122 + 2011114 + 48011123 + 36011222 + 1601133 + 1201124 

+2401223 + 80134 + 902222 + 60224 + 044 ) • (AlO) 



Combinatorial rules for generating the terms have been given before l2a|. They can be used as checks. As an example, 
we evaluate the coefficient of the term 0ii22 in -^eo- This is the number of ways of partitioning 6 objects into groups 
of 2 ones and 2 twos — 

2. Notation for traces 

Before writing out the traces we introduce the compact notation 



\n\ ■ ■p\®n2 ■ Vi® • • '1 = Tr 



(A12) 



where M'^^') is the p-th derivative of M. Also, [1 ■ p] is written as [p~\. The 'addition', ©, is not commutative, 
but all cyclic permutations of terms inside the brackets, [•••], are allowed. 'Multiplication', (denoted by the dot) is 
distributive over addition, subject to restrictions due to non-commutativity. Thus [n ■ p(Bm ■ p~\ — \_{n + m) • p~\ , but 
no simplification is possible for [n- p(Bm- p' (Bl ■ p® • ■ ■] . Traces can be added, e.g. , a [n • p] + 6 [n • p] = {a + b)[n- p~\. 
Derivatives are easy to write — 

[n-pY = -n[l®n - p] + n[{n - 1) -p® (p+ 1)]. (A13) 

The operation of taking derivatives is linear over the 'addition' in [ni ■ pi(Bn2 ■ P2® • • '1 , which is just the rule for 
taking derivatives of products. 

3. Operators 

We have the lowest orders 

01 = [11, 02 = -[2 -11 + [21. (A14) 

Then, the remaining known ones are obtained simply by applying the rules again. Since [2 • 1] ' — —2 [3 -1] + 2[1 © 2] 
and [2] ' = - [1 ® 2] + [3] , we first obtain 

03 = 2[3-ll -3[1©2] + [31. (A15) 
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Beyond the second order the results depend on the prescription for putting chemical potential on the lattice. In 
the Hasenfratz-Karsch (HK) prescription since all the odd derivatives are equal to each other, and so are the even 
derivatives, we can rewrite the above as 

03 = 2[3- 1] -3[1®2] + [1] (HK). (A16) 



a. ^th order 

At the 4th order we have 

L3-l]' = -3[4-l]+3L2-l®2], 

[102]' = -2[2- 102] + [2- 2] + Ll©3], 

[3]' = -Ll®3] + [4], (A17) 

using the rules of derivatives. Note that the coefficients in each line sum up to zero. This is a consequence of the rule 
for derivatives in eq. (|A13|) . Also note that each operator, [• • • ®ni ■ pi® • • •] , which contributes to 0„ must satisfy 
the constraint ^ riiPi = n. 

The expressions in eq. give 

04 = -6L4 • 1] + 12[2 • 1 © 2] - 3[2 • 2] - 4[1 ® 3] + [4] = -6L4 • 1] + 12[2 • 1 ® 2] - 3L2 • 2] - 4[2 • 1] + [2] , (A18) 

where the second expression holds only in the HK prescription. Note that a further consequence of the rule for 
derivatives is that the sum of the coefficients is zero for each 0„ for n > 2. 



6. 5th order 

At the 5th order we use 

[4- 1]' = -4[5- 1] +4[3-l®2], 

[2-l©2]' = -3L3-l®2] +2[le2-2] + [2- 1©3], 

[2-2]' = -2[l®2-2] +2[2®3], 

[1®3]' = -2[2- 1®3] + [2®3] + Ll®4], 

L4]' = -U©4] + L5], (A19) 
to get the following expression in the HK scheme 

05 = 24[5- 1] -60[3- 1®2] + 30[1®2- 2] +20[3- 1] - 15[1®2] + [1]. (A20) 



c. 6th order 

At the 6th order we use 

[5-l]' = -5[6-l] +5L4-1®2], 

[3- 1®2]' = -4[4- 1®2] +2[2- l®2-2] + [1®2®1®2] + [3- 1®3], 
[l®2-2]' = -2[2- l®2-2] - [1®2®1®2] + [3-2] + [1®2®3] + Ll®3®2], 
[2- 1®3]' = -3L3- 1®3] + [1®2®3] + [1®3®2] + [2- 1®4], 
[2®3]' = - Ll®2®3] + [2-3] - [1®3®2] + [2®4], 
[l©4]'==-2[2-l®4] + L2®4] + L1®5], 

[5]' = -Ll®5] + [6], (A21) 

where we see the consequences of non-commutativity of 'addition' for the first time. In the HK scheme this gives the 
result 

06 = -120[6-1] -120[4-1] +360[4-l®2] -16[2-1] +150[2-1®2] 

-180[2-l®2-2] -90[1®2®1®2] +30[3-2] -15[2-2] + [2]. (A22) 
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An earlier error in the derivative [1 © 2 • 2]' gave errors in the coefficients of [2 • 1 ® 2 • 2] (—210 instead of —180) and 
[1®2® 1® 2] (—60 instead of —90). These erroneous expressions were used in Is^. However, correcting this error 
has little numerical consequence. 

d. 7th order 

At the 7th order we need — 

[6-l]' = -6[7-ll +6L5-1®2], 

[4- 1®2]' = -5L5- 1®2] +2[3- 1®2- 2] +2[1®2®2- 1®2] + [4- 1®31, 

[2- 1®2- 2]' -3[3- 1®2- 2] - [1 ® 2 ® 2 • 1 ® 2] + 2[1 ® 3 • 2] + [2 • 1 ® 2 ® 3] + [2 ■ 1 ® 3 ® 2] , 
[1®2®1®2]' = -4[1®2®2 • 1®2] +2[1®2®1®3] +2[l®3-2], 
[3 • 1®3]' = -4[4- 1®3] + [2 • 1®2®3] + [2 • 1®3®2] + [1®2® 1®3] + [3 • 1®4], 
[3-2]' = -3[l®3-2] +3L2-2®31, 

[1®2®3]' = -2[2- 1®2®3] - [1®2®1®3] + [2-2®3] + [l®2-3] + [1®2®4], 
[1®3®2]' = -2[2- 1®3®2] - [1®2®1®3] + L2-2®3] + [l®2-3] + Ll©4®2], 
L2-l®4]' = -3[3- 1®4] + [1®2®4] + [1®4®2] + [2- 1®51, 
[2-31' = -2[1®2-31 +2[3®4], 

[2®4]' = -[1®2®4] - [1®4®2] + [3®4] + L2®5], 
[1®5]' = -2[2- 1®51 + L2®51 + L1®61, 

[61' = -L1®61 + [71. (A23) 

In the HK scheme this yields 

07 = 720[7-ll - 2520[5-l®2l + 1260[3 • 1 ® 2 • 21 - 1470[3-1®21 + 1260[1 ® 2 ® 2 • 1 ® 21 +840[5-ll 

-630[1®3-21 +420[1®2-21 +182[3-ll -63[1®21 + [ll. (A24) 

e. 8th order 



At the 8th order we need — 
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la 


)3€ 


D21' = -3[3- 1®3®21 - [2- 


1® 


3®l®2l + [l®2®3®2l + [l®3e 


D 2 • 21 + 



+ [2- 1®2-31 + [2- 1®4®21, 
[l®2®l®3l' = -2[2-l®2®l®3l-2[2-l®3®l®2l + [l®3®l®3l + [l®3®2-2l + 

+ [1®2 • 2®31 + [1®2® 1®41, 
[3- 1®41' = -4[4- 1®41 + [2- 1®2®41 + [2- 1®4®21 + [1®2®1®41 + [3- 1®51, 
[2 • 2®31' = - [1®2 • 2®31 + [1®2®3®21 + [1®3®2 • 21 +2[2®2 • 31 + [2 • 2®41, 
[1®2-31' = -3[2- 1®2-31 + [2®2-3l + [1®4®31 + [1®3®41, 
[1®2®41' = -2[2- 1®2®41 - [1®2® 1®41 + [2 • 2®41 + [1®3®41 + [1®2®51, 
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[1©4©2]' = -2[2- 1©4®2] - [1®2©1®4] + [2- 2®4] + [1©4©3] + [1®5®2], 
[2- 1®51' = -3L3- 1®51 + L1®2®51 + Ll®5®2] + [2- 1®61, 
L3®4]' = -[1®3®4] - [1®4®3] + [2-4] + [3®5l, 
[2®5]' = - [1®2®5] - [1®5®2] + [3®5] + [2®6], 
L1®61' = -2[2- 1®61 + [2®61 + L1©71, 

[71' = -L1®71 + [81. (A25) 

Putting this together yields, in the HK prescription, the expression 

08 = -5040[8-ll + 20160[6- 1®21 - 10080[4 • 1 ® 2 • 21 - 10080[3 • 1 ® 2 ® 1 ® 21 - 6720[6-ll 

-5040[2- 1®2®2 • 1®21 + 5040[2 • 1 ® 3 • 21 + 5040[1 ® 2 • 2 ® 1 ® 21 + 15120[4 • 1 ® 21 
-630[4-2l - 5040[2-l®2-2l - 2520[1 ® 2 ® 1 ® 21 -2016[4-ll +1512[2-1®21 

+420 [3 • 21 - 63 [2 • 21 - 64 [2 • ll + [21 (A26) 
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